package com.asa.numanaly;

import com.asa.HanShu;
import com.asa.HanShu2;
import com.asa.HanShu3;

/**
 * 积分计算
 * @author Administrator
 *
 */
public class D {

	
	/**
	 * 复合Simpson法则 *
	 * p175
	 * 
	 * @param a
	 * @param b
	 * @param n
	 * @return
	 */
	public static double jifen1(double a,double b,int n,HanShu hanshu) {
		
		
		
		double h = (a-b)/n;
		double xi0 = hanshu.hanshu(a)+hanshu.hanshu(b);
		
		double xi1 = 0;
		double xi2 = 0;
		
		for (int i = 0; i < n-1; i++) {
			double x = a+i*h;
			if (i%2==0) {
				xi2 = xi2 +hanshu.hanshu(x);
				
			}else {
				xi1 = xi1+hanshu.hanshu(x);
			}
		}
		
		double xi = h*(xi0+2*xi2+4*xi1)/3;
//		System.out.println(xi);
		
		return xi;

		
	}
	
	/**
	 * Romberg积分 *
	 * p185
	 * 
	 * @param a
	 * @param b
	 * @param n
	 * @return
	 */
	public static Double[][] jifen2(double a,double b,int n,HanShu hanshu) {
		
		Double[][] R = new Double[3][n+1];
		Double[][] result  = new Double[n+1][n+1];
		double h = a-b;
		R[1][1] = h*(hanshu.hanshu(a)+hanshu.hanshu(b))/2;
		result[1][1] = R[1][1];
		
//		System.out.println(R[1][1]);
		
		for (int i = 2; i <= n; i++) {
			double asa =0;
			for (int k = 1; k <= Math.pow(2, i-2); k++) {
				asa += hanshu.hanshu(a+(k-0.5)*h);
			}
			
			R[2][1] = (R[1][1]+h*asa)/2;
			
			
			for (int j = 2; j <= i; j++) {
				R[2][j] = R[2][j-1] +(R[2][j-1]-R[1][j-1])/(Math.pow(4, j-1)-1);//外推
			}
			
//			result[i] = R[2];
//			for (int j = 1; j <= n; j++) {
//				System.out.print(R[2][j]+"  "+i+"   ");
//			}
//			System.out.println();
			h = h/2;
			for (int j = 1; j <= i; j++) {
				R[1][j] = R[2][j];
				result[i][j] =  R[2][j];

			}
		}

		return result;
		
	}
	
	
	
	
	/**
	 * 自适应求积分 *
	 * 这个可以自己调整自己的步长，来处理
	 * p191
	 * 
	 * @param a
	 * @param b
	 * @param n 层次数的界
	 * @param TOL 精度要求
	 * @return
	 */
	public static Double jifen3(double a,double b,int n,double TOL,HanShu hanshu) {
		
		double[] TOLi = new double[n];
		double[] ai = new double[n];
		double[] hi = new double[n];
		double[] FAi = new double[n];
		double[] FCi = new double[n];
		double[] FBi = new double[n];
		double[] Si = new double[n];
		double[] Li = new double[n];

		double APP = 0;
		int i = 1;
		TOLi[i] = 10*TOL;
		ai[i] = a;
		hi[i] = (b-a)/2;
		FAi[i] = hanshu.hanshu(a);
		FCi[i] = hanshu.hanshu(a+hi[i]);
		FBi[i] = hanshu.hanshu(b);
		Si[i] = hi[i]*(FAi[i]+4*FCi[i]+FBi[i])/3;
		Li[i] = 1;
		
		
		
		while(i>0) {
			
			double FD = hanshu.hanshu(ai[i]+hi[i]/2);
			double FE = hanshu.hanshu(ai[i]+hi[i]/2);
			double S1 = hi[i]*(FAi[i]+4*FD+FCi[i])/6;
			double S2 = hi[i]*(FCi[i]+4*FD+FBi[i])/6;
			double v1 = ai[i];
			double v2 = FAi[i];
			double v3 = FCi[i];
			double v4 = FBi[i];
			double v5 = hi[i];
			double v6 = TOLi[i];
			double v7 = Si[i];
			double v8 = Li[i];

			
			i = i-1;
			
			
			if (Math.abs(Si[1]+Si[2]-v7)<v6) {
				
				APP+= Si[1]+Si[2];
				
			}else {
				if (v8>n) {
					return null;//算不出来了
				}else {
					
					i++;
					ai[i] = v1+v5;
					FAi[i] = v3;
					FCi[i] = FE;
					FBi[i] = v4;
					hi[i] = v5/2;
					TOLi[i] = v6/2;
					Si[i] = S2;
					Li[i] = v8+1;
					
					
					i++;
					ai[i] =v1;
					FAi[i] = v2;
					FCi[i] = FD;
					FBi[i] = v3;
					hi[i] = hi[i-1];
					TOLi[i] = TOLi[i-1];
					Si[i] = S1;
					Li[i] = Li[i-1];

					
				}
				
				
			}
			
			
		}
			
		return APP;
		
		
	}
	
	/**
	 * Simpson二重积分 *
	 * @param a
	 * @param b
	 * @param n
	 * @param hanshu
	 * @return
	 */
	public static Double jifen4(double a,double b,int n,int m,HanShu2 hanshu,HanShu dx,HanShu cx) {
		
		double h = (b-a)/n;
		
		double j1 = 0;
		double j2 = 0;
		double j3 = 0;

		for (int i = 0; i <= n; i++) {
			
			double x = a+i*h;
			double HX = (dx.hanshu(x)-cx.hanshu(x))/m;
			double K1 = hanshu.hanshu(x, cx.hanshu(x))+hanshu.hanshu(x, dx.hanshu(x));
			double K2 = 0;
			double K3 = 0;

			
			for (int j = 1; j < m-1; j++) {
				double y = cx.hanshu(x)+j*HX;
				double Q = hanshu.hanshu(x, y);
				if (j%2==0) {
					K2 = K2+Q;
				}else {
					K3 = K3+Q;
				}
			}
			
			double L = (K1+2*K2+4*K3)*HX/3;
			
			if (i==0||i==n) {
				j1 = j1+L;
			}else if (i%2==0) {
				j2 = j2+L;
			}else {
				j3 = j3+L;
			}
			
			
			
		}
		
		
		double J = h*(j1+2*j2+4*j3)/3;
		
		
		return J;
		
		
	}
	
	
	
	
	
	
	/**
	 * Gauss二重积分 *
	 * @param a
	 * @param b
	 * @param n 暂时大于1，别超过5
	 * @param m 暂时大于1，别超过5
	 * @param hanshu
	 * @param dx
	 * @param cx
	 * @return
	 */
	public static Double jifen5(double a,double b,int n,int m,HanShu2 hanshu,HanShu dx,HanShu cx) {
		
		double h1 = (b-a)/2;
		double h2 = (b+a)/2;
		double J = 0;
		
		
		for (int i = 1; i <= m; i++) {
			
			double JX = 0;
			double x = h1*rni[m][i]+h2;
			double d1 = dx.hanshu(x);
			double c1 = cx.hanshu(x);
			double k1 = (d1-c1)/2;
			double k2 = (d1+c1)/2;
			
			for (int j = 1; j <=n; j++) {
				double y = k1*rni[n][j] + k2;
				double Q = hanshu.hanshu(x, y);
				JX = JX + cni[n][j]*Q;
				
				
			}
			J = J + cni[m][i]*k1*JX;
			
			
			
		}
		J = h1*J;
		
		
		
		return J;
		
	}

	
	
	
	/**
	 * Gauss三重积分 *
	 * @param a
	 * @param b
	 * @param n 暂时大于1，别超过5
	 * @param m 暂时大于1，别超过5
	 * @param hanshu
	 * @param dx
	 * @param cx
	 * @return
	 */
	public static Double jifen6(double a,double b,int n,int m,int p,HanShu3 hanshu,HanShu dx,HanShu cx,HanShu2 ex,HanShu2 fx) {
		
		double h1 = (b-a)/2;
		double h2 = (b+a)/2;
		double J = 0;
		
		
		for (int i = 1; i <= m; i++) {
			
			double JX = 0;
			double x = h1*rni[m][i]+h2;
			double d1 = dx.hanshu(x);
			double c1 = cx.hanshu(x);
			double k1 = (d1-c1)/2;
			double k2 = (d1+c1)/2;
			
			
			
			for (int j = 0; j <= n; j++) {
				double JY = 0;
				double y = k1*rni[n][j]+k2;
				double fx1 = fx.hanshu(x, y);
				double ex1 = fx.hanshu(x, y);

				double l1 = (fx1-ex1)/2;
				double l2 = (fx1+ex1)/2;

				
				for (int k = 1; k <=p; k++) {
					double z = l1*rni[n][j] + l2;
					double Q = hanshu.hanshu(x, y,z);
					JY = JY + cni[p][k]*Q;
					
					
				}
				JX = JX + cni[n][j]*l1*JY;
			}
			
			
			J = J + cni[m][i]*k1*JX;
			
			
			
		}
		J = h1*J;
		
		
		
		return J;
		
	}
	
	
	
	
	
	
	
	
	
	
	
	
	
	//高斯求积所用到的表，这些参数如何产生的，怎么求出来的，我就不知道了，好像可以参考 参考文献StS
	public static double[][] rni = new double[6][6];
	public static double[][] cni = new double[6][6];
	static {
		
		rni[2][1] = 0.5773502692;
		rni[2][2] = -0.5773502692;

		rni[3][1] = 0.7745966692;
		rni[3][2] = 0.0;
		rni[3][3] = -0.7745966692;
		
		rni[4][1] = 0.8611363116;
		rni[4][2] = 0.3399810436;
		rni[4][3] = -0.3399810436;
		rni[4][4] = -0.8611363116;

		rni[5][1] = 0.9061798459;
		rni[5][2] = 0.5384693101;
		rni[5][3] = 0.0;
		rni[5][4] = -0.5384693101;
		rni[5][5] = -0.9061798459;
		
		
		cni[2][1] = 1.;
		cni[2][2] = 1.;

		cni[3][1] = 0.5555555556;
		cni[3][2] = 0.8888888889;
		cni[3][3] = 0.5555555556;
		
		cni[4][1] = 0.3478548451;
		cni[4][2] = 0.6521451549;
		cni[4][3] = 0.6521451549;
		cni[4][4] = 0.3478548451;

		cni[5][1] = 0.2369268850;
		cni[5][2] = 0.4786286705;
		cni[5][3] = 0.5688888889;
		cni[5][4] = 0.4786286705;
		cni[5][5] = 0.2369268850;


		
	}
	
	
	
	
	
	
	
	
	
	
	public static void main(String[] args) {
		
		HanShu hanshu = new HanShu() {
			
			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				return 0;
			}
			
			@Override
			public double hanshu(double x0) {
				// TODO Auto-generated method stub
				return Math.sin(x0);

			}
		};
		Double[][] jifen2 = jifen2(0, 3.1415926, 6, hanshu);
		for (int i = 0; i < jifen2.length; i++) {
			for (int j = 0; j < jifen2[i].length; j++) {
				System.out.print(jifen2[i][j]+"  ");
			}
			System.out.println();
		}
		
		
		
		
		
	}
	
	
	
	
}
